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ABSTRACT 

We present a new time-dependent multi-zone radiative transfer code and its first application to 
study the SSC emission of the blazar Mrk 421. The code couples Fokker-Planck and Monte 
Carlo methods, in a 2 dimensional (cylindrical) geometry. For the first time all the light travel 
time effects (LCTE) are fully considered, along with a proper, full, self-consistent, treatment 
of Compton cooling, which depends on them. We study a set of simple scenarios where the 
variability is produced by injection of relativistic electrons as a 'shock front' crosses the emis- 
sion region. We consider emission from two components, with the second component either 
being pre-existing and co-spatial and participating in the evolution of the active region (back- 
ground), or being spatially separated and independent, only diluting the observed variability 
(foreground). Temporal and spectral results of the simulation are compared to the multiwave- 
length observations of Mrk 42 1 in March 200 1 . We find parameters that can adequately fit the 
observed SEDs and multiwavelength light curves and correlations. There remain however a 
few open issues, most notably: i). The simulated data show a systematic soft intra-band X-ray 
lag. ii). The quadratic correlation between the TeV gamma-ray and X-ray flux during the de- 
cay of the flare has not been reproduced. These features turned out to be among those more 
affected by the spatial extent and geometry of the source, i.e. LCTE. The difficulty of produc- 
ing hard X-ray lags is exacerbated by a bias towards soft lags caused by the combination of 
energy dependent radiative cooling time-scales and LCTE. About the second emission com- 
ponent, our results strongly favor the scenario where it is co-spatial and it participates in the 
flare evolution, suggesting that different phases of activity may occur in the same region. The 
cases presented in this paper represent only an initial study, and despite their limited scope 
they make a strong case for the need of true time-dependent and multi-zone modeling. 
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1 INTRODUCTION 

Blazars are the most extreme (known) class of AGN. They are core- 
dominated, flat-radio-spectrum radio-loud AGN. Their properties 
are i nterpreted in terms of ra diation from relativistic jets pointing 
at us jUrrv & Padovanilll995 f). Because of relativistic beaming, jets 
greatly outshine their host galaxies thus making blazars unique lab- 
oratories for exploring jet structure, physics and origin. 

Blazars emit strongly from radio through 7-ray energies. Their 
spectral energy distribution (SE P) comprises two major contin- 
uum, non-thermal, components jUlrich et al.l 1 19971 : iFossati et alj 
[r998): the first, peaking in the IR-optical-X-ray range, is unam- 
biguously identified as synchrotron radiation of ultrarelativistic 
electrons. The nature of the second component, sometimes ex- 
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- galaxies: jets — BL Lacertae objects: individual (Mrk 421) 

- X-rays: individual (Mrk 42 1 ) — radiation mechanisms: non- 



tending to TeV energies, is less clear and under debate. It is gen- 
erally modeled as inverse Compton (IC) scattering by the same 
electrons that produce the synchrotron emission. The seed pho- 
tons can be synchrotron photons (synchrotron self-Compton, SSC 
Maraschi et al. 1992; Marscher & Travis 1996) or external radia- 
tion fields such as emission directly from the accretion disk, or 
the broad emission region (BLR), or a putat ive torus p resent on 
a larger sca le (ext ernal Compton, EC, e.g. iDermer et al. 199^; 
Sikora et al.' 'l994'; 'Ghisellini & Madai^ 1 19961 : iBlazeiowski et all 
2000 : Sikora et al. 2009). These models are generally referred to 
as leptonic models because the particles and processes responsible 
for the emitted radiation are only electrons and positrons. 



A second class of models, hadronic models, consider the role 
played by protons either by producing very high energy radiation 
directly via the synchrotron mechanism, or by initiating a particle 
cascade leading to a second leptonic population emitting a higher 
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energy synchrotron comp onent i Mannheinil \99'\ . RachenI 
'Sikora & Madejski' 2001; Arbei ter et alj |2005| ; iLevinsonI 



200C 



2ooe 



Bottcher 2007; Bottcher et al. 200^ 

The frequency of the synchrotron peak {uF^,) has emerged 
as (one of) the mo st important observ ational distinction across the 
blazar family (e.g. lFossati et alj[l998l) . leading to the classification 
of blazars as 'red' or 'bl ue' according to their SED 'color', i.e. 
the location of the pealQ. iFossati et al.l ( Il998) showed that blazar 
SEDs seem to change systematically with luminosity; the most 
powerful objects are red, while blue SEDs are associated with rel- 
atively weak sources, a result supported by s tudies of high redsh ift 
blazars and of low power B LLac objects (see Fabian et al .12001 al ibi: 



ICostamante et aU2001. ; Ghisellini et alj|201olv 

Another fundamental distinction among blazars concerns 
their 'thermal' spectral properties, where they encompass a wide 
range of phenomenology, ranging from objects with feature- 
less optical spectra (BLLac objects) to objects with quasar-like 
(broad) emission fine sp ectra (Flat Spectrum Radio Quasars, FSRQ 
IUrrv&Padovanilll995l for a review). This distinction is likely to 
have an impact on the mechanisms of production of the 7-ray com- 
ponent in different types of blazars, with BL Lacs being consistent 
with pure SSC, and FSRQ with EC. In fact in most FSRQs the 
prominence of the thermal components with respect to the syn- 
chrotron emission suggests that EC must be dominant over SSC. 
The case for BL Lacs is more ambiguous. Since particles in the ac- 
tive region (blob) in the jet would see the external radiation greatly 
amplified by relativistic aberration with respect to what we mea- 
sure, the fact that we do not directly detect any thermal component 
may not necessarily mean that in the jet rest frame its intensity is 
not competitive with the internally produced synchrotron radiation. 

However, the broad band emission of TeV detected BL Lac 
objects, like Mrk 421, is well modeled with pure SSC and stringent 
upper limits can be set on the contribution of EC to their SEDs 
dCrhisellini et al.lll998l,l2010l) . 



1.1 Variability 

Rapid and large-amplitude variability is a defining observational 
characteristic of blazars. It occurs over a wide range of time-scales 
and across the whole electromagnetic spectrum lUMch et al. 1997). 
Flux variability is often accompanied by spectral changes, typically 
more notable at energies around/above the peak of each SED com- 
ponent. Multiwavelength correlated variability studies have been a 
major component of investigation of blazars, but because of obser- 
vational limitations so far it has been focused on blue blazars. 

Blue blazars / HBLs indeed constitute a particularly interest- 
ing subclass, for their synchrotron emission peaks right in the X- 
ray band, and the high energy component reaches up to TeV 7- 
rays. The X-ray/TeV combination has been accessible observation- 
ally thanks to ground based Cherenkov telescopes and the avail- 
ability of several X-ray observatories. Hence, the brightest HBLs 
have been studied extensively. Simultaneous X-ray/7-ray observa- 
tions showed that variations around the two peaks are well corre- 
lated, providing us with diagnostics on the physical conditions and 
processes in the emission region for HBLs. 

Different models have been shown to successfully reproduce 
time-averaged or snapshot spectral energy distributions of blazars. 



^ Blue and red SED objects are also called HBL and LBL, for High or Low 
peak. 



So far, however, there has been remarkably little work taking ad- 
vantage of the information encoded in the observed time evo- 
lution of the SEDs by modeling it directly, despite the tremen- 
dous growth and improvements on the observational side, allow- 
ing in many cases to resolve SED on physically relevant time- 
scales, fueled by several successful multiwavelength campaigns 



Sambruna et al.ll2000l; iRavasio et al. 20021; Takahashi et alj 2000l; 


Fossati et al.l 


2000a b; Krawczvnski et alJ 2004 


; iBlazejowski et alj 


2005';'Rebillot et al. 20061; iGiebels et al. 20071; 


Fossati et alj 20081; 


Aharonian et al. 2009). 



The spectral time evolution has been studied and character- 
ized by means of intra- and inter-band time lags, intensity corre- 
lation, and hysteresis patterns in brightness-spectral shape space. 
The main observed features unveiled by this type of analyses seem 
to be well accounted for by attributing the 7-ray emission to SSC 
(in a one-zone homogeneous blob model), and they emerge from 
the combination of acceleration and cooling and depend on the 
relative duration of the related time-sc ales (e.g. iTakahashi et alj 
Il996l ; lUlrich et al.lll99l lKataokall2000l) . A less empirical, more 
directly theoretical interpretation of this wealth of data, requir- 
ing/exploiting the physical connection between series of spectra, 
has remained relatively basic despite the clear rich ness of the ob- 



serve d phenomenology (e.g.. iMastichiad is & Kir^ 19971; 



Dermei 



'1998'; Li & Kusunose 2000;' Sikora et alJl2001l;IKrawczvnski et al, 
i2002 ; .Bottcher & Chiana.20021) . 



1.2 Light Crossing Time Effects 

One of the biggest challenges and limitations of the current models 
comes from the dealing with light crossing time effects (LCTE), 
usually treated in a simplified way, such as simply by introducing a 
photon escape parameter (e.g., Bottcher & Chiang 2002). 

The observed variability on time-scales of hours indicates that 
light crossing time effects within the active region are very impor- 
tant and must be dealt with. There are two main aspects related 
to photon travel times that are important for an accurate study. 
The fi rst, which we can call 'external' (following Katarzvhs ki et alj 
l2008l) . is a purely geometric effect that pertains to the impact of the 
finite size of the active region on the observed variability, namely 
the delayed arrival time of the emission from different parts of the 
blob, and consequent smearing of the intrinsic variability charac- 
teristics ( Prothe roe 2002) . It is relatively simple to implement. 

The second effect, internal, pertains to the impact of these 
same delayed times on the actual physical evolution of the variabil- 
ity (as opposed to just our 'perception' of it) due to the changing 
conditions inside the active region. This effect constitutes the real 
challenge for proper multi-zone modeling. In this respect the most 
important issue is that of the photon diffusion across the blob on the 
electrons' inverse Compton losses. Proper accounting of this effect 
is significantly more complex and computationally expensive, and 
traditionally neglected under the assumption that electron cooling 
is dominated by synchrotron losses. This is, however, a strong as- 
sumption, rarely valid, as suggested by the observation that the lu- 
minosity of the synchrotron components is at best of the same order 
as the IC component, more commonly lower 

It has long been realized that a simple one-zone homoge- 
neous model is not adequate to describe the temporal evolution 
of the blazar jets, an d that LCTE must be taken into account. 
iMcHardv et~ai1 ( l2007h suggested that the observed delay between 
X-ray and infrared variations in 3 C 273 could be related to the time 
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necessary for the soft (synchrotron) photon energy density to build 
up as the they travel across the active region. 



1.3 Relevant Previous Work 

Some progress has been made to develop multi-zone models, 
though with limited success because the traditional analytical ap- 
proach requires significant assumptions, such as simple geome- 
tries or assumptions about the relevance of different physical pro- 
cesses. The inclusion of just the external LCTE is enough to yield 
new insights on SSC light curves, such as on the way the inter- 
play between cooling/acceleration time-scales and source size af- 
fects the observed light curves as a function of energy and com- 
bination of the various tim e-scales iChiaberge & Ghis ellini 19991: 
iKataoka et"Z|[2000 ; K atarzvriski et alj|2008h . In all the cited cases 
the size of the active region and the duration of the injection of fresh 
particles are related through t[^j — R/c, where R is the radius of 
a sphere or the length of a cubic region. The geometry is character- 
ized by a single length-scale. This kind of models, not accounting 
for internal LCTE and non-locally-emitted radiation for IC emis- 
sion, could yield correct results for the evolution of the electron 
distribution if synchrotron losses dominate, however even in this 
case their results for the evolution of the IC component are not re- 
alistic because they ignore the contribution of seed photons from 
othe r zones. 

ISokolov et ai] ( l2004h and ISokolov & Marsched ( l2005h were 
the first to include the internal LCTE to calculate the IC spectrum, 
for both SSC and external IC models. However, they did not prop- 
erly account for it when calculating IC energy losses. Their model 
is thus accurate only when synchrotron losses are dominant. Ob- 
servationally this corresponds (approximately) to cases where the 
peak of the lower energy component of the SED (synchrotron) is 
sign ificantly brighter th an that of the second peak (IC). 

iGraff et"aLr ( l2008l) developed a model taking into account all 
the LCTEs, but specialized to an elongated 'pipe' geometry. The 
geometry of the current implementation of their code is effectively 
one-dimensional. The lack of an actual transverse dimension repre- 
sents a significant limitation when considering the LCTE, consid- 
ering because of relativistic aberration we are effectively observing 
a jet (also) from its side. 



In this paper we introduce a more general and flexible code to 
simulate blazar variability, addressing and overcoming most of the 
limitations affecting previous efforts. 

The general features and assumptions of the code are illus- 
trated in section ij2] followed by a comparison with the results of 
other codes to test its robustness in S|3] In Sj4]we present results 
of a few test cases based on the multiwavelength observations of 
Mrk 421 in March 2001. We conclude with a discussion of this 
first application and remarks about future applications and devel- 
opments in !j5] 

In order to keep the notation light, we will use primes for blob- 
frame values sparingly, mostly to distinguish photon energies, lu- 
minosity and times (E, L, t, r). We do not prime quantities that 
are usually not ambiguous because they are only referred to in the 
blob-frame, such as magnetic field strength (B), source size (R, Z), 
electron Lorentz factor (7), density (rio). Similarly, we do not use 
primes when the context is clear (for instance in the discussion of 
the Fokker-Planck equation). 




END 



Figure 1. Basic structure and work flow of the code. The Monte Carlo block 
handles photon emission/absorption processes (e.g., synchrotron, IC, pair 
annihilation, self-absorption, escape). The MC time-step is determined at 
setup and it does not change. The Fokker-Planck block handles electron pro- 
cesses (e.g., injection, cooling, pair production, escape). The F-P time-step 
is adjusted at each iteration according to the current physical conditions. 



2 THE MONTE CARLO / FOKKER-PLANCK CODE 

Our code couples Fokker-Planck (F-P) and Monte Carlo (MC) 
methods in a 2 dimensional (cylindrical) geometry. It is built 
on the MC radiati ve transfer code deve l oped by Liang, Bottcher 
and collaborators jCanfield et ai] 1 19871 : 'Bottcher & Lian3 1200 ll : 
iBottcher et all |2003|) , parallelized by [Finke (200%. The Monte 
Carlo method is ideal for multi-zone 2D/3D radiative transfer prob- 
lems. Due to its tracking of the trajectory of every photon, LCTE 
are automatically accounted for, regardless of the geometry. We 
modified the parent code significantly in several aspects, to make it 
more general applicable, in particular to the physical conditions of 
the active region in a blazar jet. 



2.1 Code Structure 

The code separates the handling of photon and electron evolu- 
tion. The electron evolution is govern ed by the F okker-Planck 
equation, as commonly done (e.g. [ Fabia n et al. 198f :ICoDpill 19921: 
CoDpi et al." 1 993': 'Kirk et al J 1 19981 : iMaking , 1 99^ Kataoka etd] 



2000: Chiaberge & Ghisellini ll 19991) . Photons are dealt with by the 
MC part of the code, which tracks photon production and evolution 
by different mechanisms, including IC scattering with the current 
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Figure 2. The geometry of the blob model. The volume is divided into dif- 
ferent zones in r and z directions, each zone with its own electron distribu- 
tion and magnetic field. We also schematically show the setup for the vari- 
ability of the simulations presented here. The hatched layer represents the 
stationary shock ( i|2.6.U . The blob, simulation volume, is moving down- 
ward and crossing the shock front. Zones that crossed the shock at earlier 
times have had some time to radiate the newly injected energy and are plot- 
ted in lighter color shades. 



electron population, and propagation. The code's basic structure 
and work flow is illustrated in Fig.[T] 

There are two main loop structures. Since the evolution of the 
electron distribution is faster than that of the photons, each MC 
cycle contains several F-P (electron) cycles. Therefore the code 
has two main time-steps: a longer MC time-step {At'yio), within 
which the F-P equation routine performs the evolution of the elec- 
tron spectrum on shorter, variable length, time-steps (Atp_p). We 
describe them in more details in the next Sections. 



2.2 Geometry 

As illustrated in Fig. |2] the code is built with 2D cylindrical ge- 
ometry, with symmetry in the azimuthal direction. The volume has 
radius R and length Z, and it is divided evenly into zones in the 
radial and vertical directions (r and z coordinates, nr, n^). In all 
tuns presented here rir = 9 and = 30. The number of zones 
sets the resolution of the simulation for what concerns spatial in- 
homogeneities in the physical properties, either as directly set up 
or because of their different evolution (e.g. radiation energy den- 
sity will always develop a radial profile, in turn inducing a radial 
profile in the electron spectra). In the scheme adopted for this work 
the number of zones is also related to the duration of the Monte 
Carlo time-step (see §|23)- For scenarios where the variability is 
produced by a perturbation crossing the simulation volume mov- 
ing along the z axis, spatial/temporal resolution in the z direction 
is more important, hence we select a larger . For the simulations 
presented here our choice yields a Monte Carlo time resolution of 
/S.z/c/5 ^ 500 s in the observer's frame, adequate to model the 
2001 Mrk 421 data. It is possible to stud y faster variability, e.g. < 5 
minutes as observed in PKS 2155-304 jAharonian et al.l2009h . by 
increasing the number of zones at the expense of increased com- 
putational time. It is worth noting that as long as the choice of Tir, 
riz ensures that each zone is small enough to sample properly the 



variations on the physically relevant time-scales, the results of the 
simulations are substantially insensitive to these parameters. 

This geometrical setup is adequate for the case we want to 
study since the assumption is that the active region is a slice of a 
collimated jet. In principle, the code setup is flexible enough that 
slightly different geometries could be simulated via the parameter 
settings in each zone. 

Each zone has its own electron population (spectrum, density) 
and magnetic field. They can be setup individually and their time 
evolution is independent from each other, except for the effect of 
mutual illumination. Photons move freely among different zones, 
but the code assumes that electrons stay in their given zone and do 
not travel across zones; the electron Larmor radius is very small 
for the energies and (tangled) magnetic field strengths typical for 
the active region of a blazar jet, at least those studied here. For 
instance, for 7 = 10*^ and B = 0.03 G, rL = 5.7 x 10^" cm, to be 
compared with source size of the order of 10^^ cm. The radiation 
emitted by the blob is registered in the form of a pseudo -photon list 
(see ij2.3.1| >. with time, direction and energy (see also lStem et alj 

All the calculations are done in the blob rest frame. The trans- 
formation of all the quantities into the observer's frame is per- 
formed afterwards. The output is analyzed using a separate post 
processing code to produce SEDs and light curves. Since the prod- 
uct of the code is effectively a photon list, we have significant free- 
dom in the choices of bin sizes for time, energy, and angle, mostly 
limited by statistics, much in the same way as for actual observa- 
tions. Hence, we can tailor the simulation results to the characteris- 
tic of the observations that we want to reproduce (e.g. time binning, 
energy bands). 

In all cases presented here, the observed spectrum is obtained 
by integrating the beamed photons over a small solid angle cen- 
tered around the angle 6 between jet axis and observer, assumed 
as customarily to be ^ = l/F, for which also 5 = F. The typical 
width of the integration solid angle is A cos(S) ~ few x 10"^. 

2.3 The Monte Carlo section 

The MC part of the code uses the current electron distribution, as 
updated in the F-P section of the code. It includes all processes that 
involves changes in the radiation field, such as Compton scattering 
and the production of new photons by various radiative processes, 
the most important of which for our case is synchrotron emission. 
Other notable processes are pair production and annihilation, and 
synchrotron self-absorption. 

The MC time-step is currently a user-set parameter, part of a 
run input setup. It is adjusted depending on the geometry of the 
problem, e.g. shorter than the light crossing time of the smallest 
zone, and requirements of physical accuracy, for instance with re- 
spect to the fact that during each MC time-step the code does not 
change the electron distribution, which is evolved only during the 
F-P section of the code (i.e. ensuring that AtJ^j^ < t'^ooAi) for the 
highest energy electrons.) 

2.3.1 Monte Carlo particles 

Since it is impossible to follow every individual photon a common 
technique used in radiative tran sfer problems is to gro u p them into 
packe ts, pseudo-photons (e.g. lAbbott & Luc^ fl98^ : Istem eTd] 
1 19951) . to which we will refer as Monte Carlo particles. Every MC 
particle k represents rik photons with the same energy, the same 
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velocity vector, at the same position and time, carrying a total en- 
ergy of iSk = ?Jk(i^k) hvk- The nk is also referred to as statistical 
weight of the MC particle. 

The MC particles are born in the volume through emission 
processes, primarily synchrotron radiation in our case. The lumi- 
nosity of the newly radiated synchrotron contribution is computed 
and converted into MC particles with a distribution according to the 
probability given by their SED. 

The position within a given zone, time within the current time- 
step, and travel direction of the MC particle when it is generated are 
drawn randomly from the appropriate probability distributions. 

At every time-step, each MC particle moves independently, 
with some probability of being IC scattered. Absorption is handled 
as a decrease of the statistical weight of the MC particles. 

When a MC particle reaches the volume boundary, it is 
recorded with the full information of the escape time, position, di- 
rection, and energy, forming a list of emitted photons. 



2.4 The Fokker-Planck equation 

In each zone, the temporal evolution of the local electron popula- 
tion is obtained by solving the Fokker-Planck equation: 



dt 



_d_ 

1 I" 



iV(7,f)7(7,f) 
Ni^,t)D{^,t) 



+ Q{'f,t)- 



(1) 



tc 



N{'y, t) is the electron spectrum, 7 is the random Lorentz factor of 
electrons, 7 is the total heating/cooling rate. The IC cooling uses 
the time dependent radiation field calculated in the Monte Carlo 
part of the code, with LCTEs accounted for, which is considered 
constant for the duration of the F-P section of the code. The full 
Klein-Nishina (K-N) scattering cross section is used (see i)2.6.3t . 
D{-y, t) is the dispersion coefficient which is not important for the 
type of scenarios presented in this work, and it then set to zero. For 
generality, the -D(7, t) term is still included in the solving of the F- 
P equation. Q(7, t) is the electron injection term. Because, as noted 
in §|2j2] the electrons' Larmor radius is much smaller than the size 
of the simulation zones, the particle escape term is not considere(|f]. 

The time-step of the F-P loop is adjusted automatically de- 
pending on the rate of change (gain or loss) of energy of the parti- 
cles to ensure a physically meaningful solution. It is constrained to 
be shorter than one fourth of the MC time-step. 

Rather than u s ing th e discretiza tion scheme p r opose d by 
iNavakshin & Meli j l ll998l) . as done in ISottcher et alj ( l2003h . we 
choose to adopt the implicit difference scheme proposed by 
IChang&Coopeil ( 1197( f). This scheme guarantees non-negative so- 
lutions, which in runs with the original scheme resulted in wild 
oscillations of the electron distribution at the high energy end 
(for a discussion of t his issue please refer to the appendix of 
IChang&Coopeilll970l) . 

The energy grid used for the electrons is logarithmic in kinetic 
energy Xj = 7^ — 1, with 200 mesh points from a^min = 0.18 to 
a;max = 3.1 X lO'', i.e. x^ 



2 Except for the test runs discussed in Section i|3.2| for consistency with the 
model with which we are comparing the results. 



After rewriting l|T) as 

dN{-i,t) _ d 



dt d-y 



It is possible to discretize it as: 



2 dj 

+ Qii,t) 



7V(7,t) 



(2) 



At 



1 

Ax, 



+ 



+ 



Axj+i/a 
1 

1 



Cj+1/2 Wj+1/2 



AXj_l/2 
1 

Ax,_i 



C3-I/2 -j- 



1 _ e^-™3-i/2 / ' 



C,_i/2 W^,_i/2 Nf+,' 



tesc 



(3) 



with 



Bj+1/2 = -T^tij +7j+i] + 



Dj+i - Dj 



Cj+1/2 = -{Dj + Dj+i) 



+ — AXj.|_i/2-Bj + l/2/C'j + i/2, 
W.j + i/2 = Wj + i/2/(exp(Wj + i/2) - !)■ 

Here the j ± 1/2 subscripts refer to quantities computed as the 
average values of the two adjacent grid points, such as 

Cj+1/2 = -^(Cj + Cj+i) 

An exception is that of 



Axj = ^Aa:j+i/2Aa;j_i/2 
with 

^Xj + i/2 = Xj+l - Xj 
'^^j-1/2 = Xj — Xj-\. 

In order to avoid infinity in our calculation, we set D — 10~*°s~^ 
instead of D = 0. 

The tridiagonal matrix fo rmed by equation ([3} can be solved 
using the standard algorithm in IPress etal.l ( [T992h . 

2.5 Synchrotron and inverse Compton 

The synchrotron spectrum is calculated adopting the single particle 
emissivi ty averaged over an iso t ropic d istribution of pitch an gles 
given bv lCrusius & Schlickeiseil l ll986h jGhisellini et al.lll988h : 
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(4) 



where ctt is the Thomson cross section and 

E eB V 



1 = 



27r rricC 



y ■ 



with E the total electron energy, and Kx (y) is the modified Bessel 
function of order x. 

The total emitted synchrotron power and self-absorption 
coefficients are calculate d according to the formulte in 
iRvbicki & Lightmad (l91^ . 

J2 



■J dEP{v,E) E'^ 



N{E - hv) N(E) 
{E-hvY 



£2 



(5) 



where P{v, E) is the synchrotron spectrum given above (O. 

For the total Compton cross section, we used th e angle- 
averaged cross section given in lCoppi & Blandforcj ( Il99(]|) : 



a(a;,7) 



Sctt 



8 



X 1 
2 ^ 2(1 + 



- + 



9 + X + - ) ln(l + x) + 'iLi2{-~x) 



x=2'^{l+P)u (6) 



Where 1/12 (z) is the dilogarithm, which is evaluated numerically. 
To get the total cross section for a photon in an electron medium 
we need to integrate over 7, weighted by the electron energy distri- 
bution. 

2.6 Other major changes 

Besides changing the numerical scheme to solve the F-P equation, 
we implemented several other major changes in the code. 



2.6.1 Injection of electrons 

The model of the electron injection process, as implemented cur- 
rently, involves a stationary shock perpendicular to the axis of the 
cylinder (jet) (Fig.O. Hence, in the frame of the blob the shock is 
traveling across the blob with a speed equal to the bulk velocity of 
the blob Ubuik ^ c. This scenari o is similar to the one discussed 
bv IChiaberee & Ghisellinil (IT999h. The thickness of the shock is 
treated as negligible, in the sense that it is considered active only in 
one zone at any given time, i.e. it never splits across a zone bound- 
ary. However, during the time it takes to travel along a Z-zone, 
Az/c, particles are injected in the entire zone volume. From this 
point of view the 'shock' thickness is not negligible. Provided that 
the of each zone is small this approximation is reasonable. As 
noted, for the cases presented here, Az/c/5 ^ 500 s. The total 
duration of injection is thus fj^j = ^/wbuik, and each slice of the 
simulation volume along the z axis will eventually have an injected 
energy of L[„j Az/iibuik, where Az = Z/n^ is the thickness of 
one slice. 

Electron injection is included in the Fokker-Planck equation 
through the term Q{-f,t). The shock moves at the speed of «buik 
every F-P time-step. When the shock front is located in a given 
zone, electron injection is active (Q 7^ 0), otherwise Q = 0. In the 
simulations presented here the injected electrons have a power law 
distribution with a high energy exponential cutoff 



Q(7) 



70 



-7/7ir 



7 > 7mm, 



The value of the normalization Qo is controlled by the parameter 
r' 



Injection and acceleration time-scales and durations are in 
principle independent from other parameters and could be set di- 
rectly on the basis of a hypothesis on the details of physical pro- 
cesses underlying them. In this work we are treating injection, and 
in turn the implied process for accelerating the newly injected par- 
ticles, phenomenologically, affording ourselves the freedom to as- 
sume their spectrum and time-scales. 

The underlying physical mechanism for the injection process 
is not specified. First order Fermi acceleration at a shock front or 
second order Fermi acc eleratio n by a plasma turbulence are two 
possible processes (e.g. Drurv 198 3: Blandford & Eichler 198 

Kirk et al.lll998l : lKatarzvnski et 



Gaissedl 199 ll : |Protheroelfl996l : 
2006h . 



2.6.2 Splitting of MC particles 

A major difficulty of using the Monte Carlo method to model 
broadband IC emission, in the physical conditions typical of blazar 
jets, is the low pseudo-photon statistics at high frequencies. Obser- 
vations are affected by a very similar problem. 

Blazar SEDs are approximately flat in {vF^) over wide range 
of energies. In blue blazars typical energies for the electrons re- 
sponsible for the (vFv) emission peaks, occurring in UV-X-ray and 
TeV bands, are 7 ~ 10" - 10^ When a photon (for us a MC par- 
ticle) is scattered to the TeV range, the energy of that MC parti- 
cle will increase by about 9-1 1 orders of magnitude depending on 
whether it was an X-ray or optical photon, and its 'flux' will de- 
crease by the same factofl making the statistics of the high energy 
IC component very poor. 

An additional challenge that we face is that the IC scattering 
probability is very small. Under most reasonable conditions the ac- 
tive blob is very optically thin. 

In order to mitigate these problems, we introduced a method 
relying on the splitting of MC particles. The basic idea is that since 
every MC particle represents a packet of real photons treated to- 
gether, it is always possible to divide them into smaller packets. If 
this splitting is applied in appropriate conditions, it is possible to 
achieve a reasonable statistics on MC particles at high frequencies 
with reasonable computer resources. 

We have implemented MC particle splitting in three different 
instances within the context of the computation of IC scattering. 

(i) The first splitting is applied to every MC particle when it is 
considered for IC scattering. It is split into a large number of iden- 
tical subparticles (e.g. ~ 10'^). The choice of this number depends 
on the trade off between improving the statistics of the high energy 
photon spectrum and cost in terms of computing resources (time 
and memory), and it was based on empirical testing. Whether a 
particle is scattered or not is determined by comparing the distance 
it would travel with a distance to the next scattering stochastically 
determined from its mean free path. Every MC subparticle draws 
a separate random number, and in turn has its own probability of 
being scattered. All non-scattered MC subparticles are recombined 
into a MC particle, and travel to a new position. The subparticles 

^ For constant statistical weight, the discretized spectram would have 
Ni ~ N(vi){l\v)i MC particles in each bin. Our grid of photon en- 
ergies is equally spaced logarithmically, so we can rewrite it as Ni ~ 
N(vi) Vi (A In v)i. where (A In v^i = A In is a constant. Hence for 
a photon spectrum N{u) on , the relative statistics of our discrete pho- 
ton spectrum goes like For an approximately flat SED, i.e. F ~ 2, 
this goes like !/ r ^ . 
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that do scatter (usually a small number) will be scattered separately, 
to independent energies and directions (but see below). This first 
splitting does not necessarily save a computational time, but it de- 
creases dramatically the memory allocation requirement to achieve 
the desired statistics at the highest SED energies. 

(ii) The second instance of splitting is applied to MC 
(sub)particles that are being scattered. They are divided into an- 
other large number (e.g. ~ 10^), and each of these MC sub- 
subparticles will be scattered separately, to a frequency and direc- 
tion uncorrelated with those of the other particles. This splitting 
allows us to concentrate computation cycles on the rare events of 
scatterings, which is what we are most interested in. 

(iii) Even with this second splitting, at highest energies the 
statistics of the IC photon spectrum remains very poor. To alle- 
viate this problem, we implemented a third instance of MC particle 
splitting. It is triggered when one of the already twice-split MC par- 
ticles is scattered to very high frequency, above a threshold that is 
set a priori and constant for each run, tailored to the characteris- 
tics of the studied SED. This MC particle is split again, and each 
of its subparticles is re-scattered from the original frequency. That 
scattering is accepted only when the scattered frequency is higher 
than the preset threshold, otherwise it goes back and draws another 
random number. This third splitting offers the benefit of avoiding 
the use of a much larger number of subparticles in the second in- 
stance of splitting, and subsequently avoiding the production of a 
very large number of MC particles to be recorded in the computer 
memory. 

Splitting causes the number of MC particles to grow during the 
simulation. Nevertheless, the advantage over directly setting up the 
simulation with more MC particles is significant both in terms of 
number of MC particles and more importantly because the new MC 
particles are created where they are most needed, thus increasing 
greatly the efficiency of the code. In typical runs the increase in 
the number of MC particles due to the splitting is modest, of the 
order of 10-20 per cent of the number of newly emitted synchrotron 
photons at each MC step. 

2.6.3 Arbitrary electron energy distribution 

Although the F-P equation can calculate the time dependent evo- 
lution of the electrons with arbitrary spectrum, earlier versions of 
the code forced the decomposition of the electron population into a 
low-energy thermal population plus a high-energy power law tail. 
The emissivity of cyclotron, non-thermal synchrotron and thermal 
bremsstrahlung radiation processes were calculated on the basis 
of this decomposition. The calculations of the synchrotron self- 
absorption coefficient and the total scattering cross section of a 
photon in the medium were dependent on this 'thermal plus power 
law' approximation as well. In order to make the code more gen- 
eral, and in particular more suitable for blazar simulations, in which 
there is usually a dominant non-thermal lepton population, we have 
entirely rewritten the relevant sections of the code. The code now 
calculates all physical quantities using the actual electron spectrum, 
as updated by solving the F-P equation (see m.5\ 

2.7 Deactivated Features 

Some features of the code have been deactivated in this study. 
Among these are the cyclotron and bremsstrahlung emission and 
Coulomb scattering of electrons with protons, all considered not 
important in blazar jets. Others are turned off because they are not 



the focus of this paper; these include external sources of photons, 
which will be subject of future investigations. 



3 TEST RUNS 

In order to test the reliability and robustness of the code, we com- 
pared the results of our code with those of other authors using dif- 
ferent codes, for cases where the codes' capabilities are compara- 
ble. We first compare the results with a non time-dependent code 
to test the MC radiative part of the code. Then we compare the 
electron evolution with a time dependent code, in a single-zone ho- 
mogeneous case. Generally the results match very well. 

3.1 Steady state SED of homogeneous models 

To test how the code handles the radi ative processes, we t ry to re- 
produce the theoretical SED shown in Fossati et al.l ( 1200 8). for the 
non-extreme parameter choice (solid line in their fig. 10a). That 
SED was computed with a single-zone homogeneous SSC model. 
The electrons are assumed to be continuously injected and reach a 
steady state (e.g. Ghisellini et al. 1998). For this test, we take the 
equilibrium electron distribution calculated in the homogeneous 
code as our input electron distribution, and turn off the F-P evo- 
lution of the electrons. We cut the volume into several identical 
zones just to make use of the parallel structure of the code. Since 
the single-zone model uses a spherical geometry, while our MC 
model uses a cylindrical geometry, we choose to use the same ra- 
dius (R = 10^" cm), but with the height of the cylinder Z = 4/3i?, 
in order for the two models to have the same volume. The produced 
SED is shown in Fig. [3^ in t he top panel as black histogram, di- 
rectly compared with that of iFossati et al. I l2008i) . In general the 
two SEDs match well, except for a slight discrepancy around the 
peak of the IC component. This arises from the fact that the sin- 
gle zone model uses a step function to approximate the K-N cross 
section, while our code implements the full K-N cross section. 

We then also tested our code with the step function approxi- 
mation. The result in shown in Fig. [3^, middle panel. The overall 
shapes of the SEDs match better. The total luminosity seems a little 
higher in the MC model. However, it is worth noting that although 
we are matching the volume, the geometry is different in the two 
codes and this has a small effect on the IC component. Moreover, 
in order to achieve a reasonable statistics the emitted photons are 
integrated over a finite solid angle, i.e. a range of angles 6 (photon 
direction angle with respect to the jet axis, in the observer frame). 
Hence for a given bulk Lorentz factor Pbuik = 26, we are effec- 
tively integrating over a range of Doppler factors 5 ~ 23.5 ~ 28.7, 
not exactly 5 = 26 as for the comparison model. 

3.2 Temporal evolution (one zone model) 

The other important aspect of our MC/F-P code, the Fokker- 
Planck evolution of the electrons, was tested by comparing the 
co de with the one-zone time dep endent homogeneous SSC code 
by IChiaberge & GhiselUnil ([199^. We set the number of zones 
to one, and used a power law injection with the same parame- 
ters they used (B = IG, 7min = 1, 7max = 10^, p = 1.7, 
iJnj = 3.69 X 10"^ erg/s, t^sc = 1.5-R/c, tj^j = R/c), except 
that our geometry is a cylinder with R = 1.1547 x IQ^^ cm, 
Z = IQ^^ cm, while they use a sphere with R = lO^^cm. 

The electron spectra at different times are shown in Fig. [SJj; 
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Figure 3. (a) Results of the steady state homogeneous model simulation. The smooth red lines are the SEDs from lFossati et alj i2008l) . The black histoeram is 
the SED produced by our MC/F-P model run with the closest possible source setup. In the top panel we used the full K-N cross section, while in the middle 
we used a step function approximation as done in Fossati et al. (2008). The bottom panel shows the fractional deviation, Ay/{y). between these two latter 
SEDs. Positive sign corresponds to the one-zone spectrum having a larger value, (b) Results of th e time-dependent homogeneous model simulation. Electron 
distributions, as 7^ Af (7), for different times for our MC/F-P code (top), or the one-zone code of lChiaberge & Gtilsel lini ( 1 999) (middle). In different colors 
and line types we plot the spectra at the following times (in units of R/c): 0.5 (magenta, long dash dot), 1 (black, solid), 1.25 (red, dotted), 1.5 (green, sort 
dash), 2 (blue, long dash), 3 (cyan, dot dash). In the bottom panel we show the fractional deviation. Ay/ (y), between the spectra for the two codes. Lines are 
tnincated at the energy where 7'^Af-y drops a factor of 30 below its peak value. 



the upper panel shows the one produced by the MC code, the mid- 
dle panel the one produced by the one-zone code, while the bottom 
panel shows the deviation. The two spectra match reasonably well, 
giving us confidence that our code handles the evolution of the elec- 
tron distribution correctly. 



4 APPLICATION TO MRK 421 

Mrk 421 is the archetypical 'blue' blazar, the most luminous and 
best monitored object in the UV, X-ray and TeV b ands. It was 
the fi rst extragalactic source detected at TeV energies JPunch et aD 
I1992I) . As such it has been the target of multiple multiwave- 
length campaigns with excellent simultaneous coverage by X-ray 
and Te V telescopes. (Takahashi et al. 2000; Maraschi et al. 1999; 
Fossati et al. 2000a, b; Kra wczvnski et alj|20 01; Blazeiowski et al. 
2005.; ,Rebillot et al..,2006, ; jGiebels et al.ll2007.; .Fossati et al.ii2008i ; 
Donnarumma et alj2009l) . 

For a first application of our code we focused on one of best 
flares ever observed, occurred on 2001 March 19 (Fossati et al. 
l2008l) . It was a well defined, isolated, outburst that was observed 
both in the X-ray and 7-ray bands from its onset through its 
peak and decay. It uniquely comprised several rare favorable fea- 
tures, namely absence of data gaps (except RossiXTE's short or- 
bital gaps), excellent TeV coverage by the HEGRA and Whipple 



telescopes and large amplitude variation in both X-ray and 7-ray 
bands. 

4.1 Observational constraints and goals 

We aim to reproduce several observational features. Some of them 
can be regarded as constraints on the setup of a basehne model, 
as they provide guidance on the general properties and para meter 
value s yielding an acceptable fit to the SEDs (e.g. ( Tavecchio et alj 
Il998l iBednarek & Prothe"ro3ll99l and iFossati eTal... 2008. for an 
example specific for the observations studied here). In this respect, 
we have five fundamental observables we want to match: 

• The peak frequencies of the synchrotron and IC components, 
t'p.s, i^p,ic, which for Mrk 421 are observed in the X-ray and 7-ray 
bands. 

• The peak luminosity and the relative strengths of the two SED 
components, vLp^s, vLpjc- 

• The variability timescale (tvar)- Combined with an hypothesis 
on the Doppler factor it provides a constraint on the size of the blob. 
For Mrk 421 in X-ray and 7-ray it is typically of the order of tens 
of kiloseconds. 

Besides giving an indication about the size of the active re- 
gion, the latter can be different for different energy bands and in 
turn its energy dependence can provide additional constraints on 
the model parameters and source geometry. 
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There is then a set of observational features whose explana- 
tion remains to a large extent an open question. They represent the 
ultimate goal of our work and the driver for the development of a 
time dependent multi-zone model. 

(i) The quasi-symmetry of flare light curves, showing sim- 
ilar rising and falling time-scales, both in X-rays and 7-rays. 
The symmetry seems to be a quite common feature at several 
wavelengths. It would seem to support the interpretation that the 
flare evolution is governed by the geometry of th e active region 
dChiaberge & Ghisellinil l 1999i - |Kataoka et al.ll200( j). However, this 
could be true only if all other (energy dependent) time-scales are 
shorter than the blob crossing time, or relevant geometric time- 
scale, or only for emission at energies for which this is true. 

(ii) The characteristics of the multiwavelength correlated inten- 
sity variations. The flare amplitude is generally larger in 7-rays than 
in the X-ray band, and flux variations show a quadratic (or higher 
order) relationship that holds during both the rise and the decay 
phases of the flare. This behavior was observed in Mrk 421 on 2001 
March 19, a nd also for other 'clean ' flares, including for other blue 
blazars (e.g. lAharonian et"al]|2009l . reporting a cubic variation for 
PKS 2155-304). 

(iii) The existence and length of inter-band (X-ray vs. 7-ray) 
and intra-band (soft X-ray vs. hard X-ray) time lags, often with 
changing sign from flare to flare (see references giv en above for 
Mrk 4 21). In the isolated flare of 2001 March 19, iFossati et all 
j2008l) report a possible lag of about 2 kilo-seconds of the TeV 
flux with respect to a soft X-ray band (2^ keV), whereas TeV and 
harder X-rays (9-15 keV) were consistent with no lag. In turn an 
X-ray intraband lag was detected. 

(iv) The fact that even during large outbursts the optical flux 
changes little. This may constrain the characteristics of the par- 
ticle injection, such as their spectrum (shape and density) and 
energy span. On the other hand, the time dependent spectral 
behavior of blazars has led people to speculate that there is 
more than one component contributing to the blazar emission 
iFossati et al.' 2000b; Rraw czvnski et alll2004l : iBlazeiows ki et alj 
[2005; Ushio et al. 2009). It is not clear if this additional zone is 
co-spatial with the zone undergoing the flare or it is far enough 
elsewhere along the jet that the two do not interfere with each other 
and evolve independently. 

(v) SED shape, and its time variations, particularly around the 
two peaks. For Mrk 421 we mostly focus on the X-ray and TeV 
7-ray spectra. 

These features have been observed in several instances for Mrk 421 , 
mostly cleanly in the case of the 2001 March 19 flare, and the 
other well studied TeV detected blue blazars. For the brightest blue 
blazars there is an extensive database of multiwavelength obser- 
vations and studies of time resolved spectral variability. The phe- 
nomenology is richer and more complex than the few items just 
introduced, on which we focus. In this respect, one of the most 
interesting findings is the observation of a correlation between 
luminosity and position of the peak of the synchrotron compo- 
nent (e.g.rTav ecchio et al. 2001; Fossati et al. 2000b ; Tanihat aet alj 
l2004lTramacere et al.ll2007l) . 

In this work we are mostly aiming at illustrating the capa- 
bilities of our code with respect to investigating the above obser- 
vational findings, by presenting the results of simulations of three 
simple scenarios. 



4.2 On model parameters 

This code affords us a great freedom. In particular, we could setup 
each zone with different initial conditions. However, for the sce- 
narios presented in this work we took a conservative approach and 
setup each zone with identical values for the usual set of physical 
parameters. 

Our homogeneous (at least initially) SSC model is defined by 
the following quantities (see also TablelTJ: 

• source size/geometry (7?, Z, or aspect ratio), 

• Lorentz factor (F), 

• magnetic field strength (B), 

• various parameters describing the electron spectrum, e.g. four 
for an injected power law: 7min, 7max, p, L[^j . For a broken power 
they would be six because there would be a spectral break 7b and 
two spectral indices (pi, P2) instead of one. 

With simple considerations we can reduce the number of 
model parameters to constrain from 8 (or 10) to 5 (B, F, R, 7max, 
L[^j) and as illustrated in the previous section we have 5+ funda- 
mental observables to do it. 

The source aspect ratio can be at least qualitatively constrained 
by the profile of the flare light curve, for in first approximation 
extreme geometries would yield fairly distinctive flare shapes due 
to LCTE. For this work we adopted a conservative, stocky, volume 
aspect ratio R/Z — 3/4, i.e. width:length = 3:2. 

Among the electron spectrum parameters, 7min and p (or pi) 
can be set with reasonable confidence based on considerations on 
the SED shape and variability (or lack thereof) at frequencies below 
the synchrotron peak. The precise value of 7min is however not well 
constrained by observations. The emission by electrons at 7 10'^ 
would be below the optical band, where there is not much simul- 
taneous coverage, and emission by much lower energy electrons 
would fall in a band (i.e. u ^ 10^^ Hz) where observations sug- 
gest that the SED is dominated by radiation f rom other regions of 
the jet (e.g. iKellermann & Paulinv-Tothlll98li) . Moreover, cooling 
time-scales for electrons at those energies are long compared with 
the typical duration high-density multiwavelength campaigns (see 
eg .11 lb. making it difficult to set a constraint on 7miii based on vari- 
ability. Higher values of 7min affect the synchrotron emission in the 
optical band and in turn the IC component, mainly in the GeV band, 
and therefore we can assess their viability with curre nt and future 
observations. Given that during the 2001 campaign jFossati et alj 
,2008.) there seemed to be a modest level of variability in the optical 
band. Amy — 0.4, though not directly from observations simul- 
taneous with the March 19 flare, we simulated scenarios where the 
injected electron population has a relatively low 7min = 50 (Ta- 
ble[Tll. We choose to truncate the electron distribution at this value 
also because the number of low energy electrons grows rapidly, 
thus increasing significantly the computational time without adding 
much to the investigation presented here; as noted, emission from 
lower energy electrons would not be detectable, and they would not 
significantly alter the properties of the emission and its variability 
observed in blue blazars. This is of course an assumption that is 
valid for this work but that should be revisited, for instance for the 
study of red blazars. 

The spectral indices of the injected electron distributions p or 
Pi (P2 = pi + 1 as expected for a cooling break) are mostly con- 
strained by the shape of the synchrotron SED at energies below 
the optical range. The preferred value forp,pi = 1.5 constitutes a 
somewhat hard spectrum but it is consistent with values discussed 
by several particle acceleration studies. In particular, stochastic 
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Table 1. Summary of model parameters 



Case general source parameters back-/fore-ground component injected component 

R Z T B 7„,in 7b 7max "c 7min 7max L[- 

1016 1016 1040 

cm cm G cm~^ ergs"'^ 

1: with 'background' 1.0 1.33 33 0.1 50 2 x 10'* 2 X 10^ 4.0 50 1.9 X 10^ 5.5 

2: with 'foreground' 1.0 1.33 33 0.08 50 1 X lO"* 1 X lO'^ 6.0 50 1.9 x 10'^ 6.0 

3: better TeV spectrum 1.5 2.0 46 0.035 50 2 x 10"' 2x10^ 1.56 50 1.9 x 10^ 3.2 

All background or foreground components electron spectra are broken power laws (with exponential cutoff), with spectral indices pi = 1.5, P2 = 2.5 
(= pi + 1). The injected power law has spectral index p = 1.5 in all cases. 



(2'"' -order Fermi acceleration) and acceleration at relativistic shear 
layers have been suggested to produce very hard (p < 2) parti- 
cle spectra (''Virtanen & Va inidl2005l : IStawarz & Ostrowskill2002l : 
iRieger & D uffv 2004, 2006). 

In the context of this discussion, it is worth adding that we 
considered scenarios with and without a pre-existing population of 
relativistic electrons or an external 'diluting' SED (see ij |4.4l and 
§|43]for details), and their characteristics could be regarded as an 
additional degree of freedom of our modeling. In this respect, how- 
ever, while a particular choice of values has some effect on the best 
parameters for the component responsible for the variability, its ef- 
fect is fairly limited and the most important point is the existence 
or not of such secondary component (see §|5}. 

Next, we illustrate some general arguments and estimates for 
values of the fundamental physical parameters. We then present 
and discuss the results obtained with the set of parameters that we 
deemed more successful, and in turn 'fit' the SEDs and light curves 
of the 2001 March 19 outburst testing several different parameter 
combinations, including the possible dilution by emission from a 
different region of the jet not involved in the flare, and the presence 
of a pre-existing electron population in the region that becomes 
active. 



4.3 Estimates of active region parameters from observables 

Key parameters in the modeling of blazars with the SSC model in- 
clude the Lorentz factor, the magnetic field strength, the size of the 
volume, and the energy of the electrons that are responsible for the 
synchrotron peak of the SED, 7p. This latter is associated with a 
break in the electron distribution or its maximum, depending on 
the spectral index. We use the observational results of iFossati et alj 
1I2OO8) as the benchmark of our analysis. There are several ob- 
served features that constrain the value of these parameters (see 
previous section). Additionally, independent estimates of the rela- 
tivistic beaming parameters of blazars, from observed superlumi- 
nal motion as well population statistics, yield Lorentz factors of 
the order of tens (.Urrv & Padovani 1995). As we mentioned be- 
fore, we make the common assumption to be observing the s ource 
at the angle 9 = l/V, hence 5 = T (see lCohen et al.ll2007 l. for a 
deeper statistical analysis, showing that the most likely combina- 
tion is T sin 6 ~ 0.7). 

The observed peak of the synchrotron component (at en- 
ergy i?p,s) results from the combination of electrons' 7, B and 
5. Assuming mono-energetic emission the relationship is iSp.s = 



I'B 7p 5. For Mrk 421 i5p,s — 1 keV. Parameterizecj^ on fiducial 



4 Because the redshift of Mrk 421 is small, 
left out factors (1 -I- 2). 



0.031, for simplicity we 



values for these three parameters the 7p of the electrons emitting at 
the synchrotron peak is: 



7p ~ 1.7 X 10^ 



B 



0.1 G 



-1/2 , ^ ^-1/2 ^ ^^^^ ^1/2 



30 



IkeV 
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If the IC component peak resulted from scattering of photons of 
the synchrotron peak in Thomson regime we could directly infer 
the energy of the electrons emitting at both SED peaks as 



7p 
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P,IC 
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P,S 
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1/2 



P,S 
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(8) 



with i5p,ic is the peak ene r gy of the IC component. However, as 
discussed by iFossati et alj ( 20081 . see Fig. 1 1 therein), the SED 
shape and variability time-scale observed in Mrk 421 in 2001 fa- 
vor parameters such that the scattering between electrons at 7p and 
synchrotr on peak photons at E^ s would happen in the K-N reg ime 
(see also iTavecchio et ailll99 i; iBednarek & Protherod 11993) . In 
this case the IC peak energy would be largely independent of i5p,s 
and the expression would instead be: 



7p 



P,IC 



(9) 



Requiring that the condition for Thomson regime, 7^' ^3/4 
(where x' = -E^argot/l'TioC^)), holds true for S^rgct = ^^p,s and 
7 = 7p, one can derive a rough estimate of what {B, 5) combi- 
nation would be necessary to push into the Thomson regime the 
scattering between 7p and its own synchrotron photons, emitted at 
Ev,s- 



B 



^ 220 



E, 



P,S 



(10) 



0.1 Gy V3oy VlkeV, 

As shown by IFossati et alj feoOSh . it is indeed possible to achieve 
an acceptable SED fit with high B and 5. However, while this kind 
of model matches equally well a static SED, its smaller blob size 
and extreme Lorentz contraction make it implausible when com- 
pared with more dynamic observational findings, beginning with 
the variability time-scales. 

The rest frame synchrotron cooling time can be expressed as 
function of electron energy and the magnetic field: 



7.7 x 10" 

7^2 



(11) 



or, more directly related to observables, in terms of observed pho- 
ton energies: 



Tcooi.s = 4.6x10" 



0.1 G 



30 



-1/2 
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or, in the observer's frame, 

0.1 G; Uoj VlkeV 



Tcool.S ~ 15.1 



Es 



-1/2 



ks (13) 



showing its dependence on the inverse square root of the energy of 
the observed photons. 

A general constraint among the observed variability time- 
scale and source size and time-scale of the acceleration, injection 
or cooling process is: 



It ' 

max I Tcool 1 "^acc 1 t-in 



R Z 



(14) 



As noted in Section |2.6.1| in this work we take a simplified ap- 
proach, whereby we do not specify the acceleration mechanisms 
underlying the particle injection, and we choose to link the injec- 
tion time-scale to the geometry of the source, namely its dimension 
along the jet axis, Z. Hence we have a simplified relationship with 
the observed variability, and considering that the 2001 March 19 
flare has a flux doubling and halving time of the order of 10^ s, we 
have approximately 



R Z 



c c 



R 



10'' (5 s 



(15) 



Please note that this constraint could actually vary with the ob- 
served band because some time-scales are likely to be energy de- 
pendent. 

If the IC cooling rate is similar to the synchrotron cooling rate 



' '''cooi.s/2- The condition t'^ooi < R/c translates into 



Es > 0.46 



R 



1016 cm 



B 



0.1 G 



30,' 



(16) 



From the constraints and relationships illustrated above we in- 
fer that a good starting point to model the SED of Mrk 42 1 is a com- 
bination of i? ~ lO^'^ cm, B ~ 0.1 G, r ~ 5 ~ 30, 7p ~ 10^ 

Because of computational limitations we did not perform an 
actual fit to identify the best set of parameters values reproducing 
the SED and the flare evolution properties. We explored a limited 
range of values for R, B, T around the values obtained from the 
above analysis, and focused on changes of the maximum electron 
energy 7max and the injected power Lj^j . 

We ran a large number of short simulations aimed at sampling 
a reasonable range of values around our initial guesses and evalu- 
ated them mostly on the basis of their matching the X-ray spectra 
and variability. In a second stage we honed in on the best cases and 
adjusted the parameters by means of full length simulation^ 



4.4 Case 1: injection in a blob with a pre-existing 
(bacliground) electron population 

In all cases presented in this paper, the outburst is attributed to the 
injection into our active volume of a new population of higher en- 
ergy electrons, with fixed injected spectrum (power law with expo- 
nential cutoff). 

In this first scenario the blob is not empty, but it is filled with 



A typical run takes around 24 hours on eight Xeon 2.83 GHz CPU cores, 
using up to 16 GB of memory. As currently implemented the code does 
not scale well with the number of CPUs, only gaining a factor of three in 
speed by going to 96 CPUs. The bottleneck is mostly due to the longer 
computational time required by the zones with larger volume because it 
scales with the number of photons contained in each zone. 



a 'background' population of electrons, homogeneous throughout 
the volume. These electrons serve as a slowly evolving component 
in the electron distribution and in turn the SED, which can be re- 
garded as the remnants of a previous phase of activity. They partic- 
ipate fully in the time evolution of the blob, cooling and emitting 
radiation. 

The overall best case has the following parameters: R = 
10^** cm, Z = 4/3 X lO^'^cm, B = O.IG, T = 33. Parameters 
for this and all following cases are summarized in Table[T] 

At t = the electron spectrum for the 'background' popula- 
tion is a broken power law distribution: 



iV(7) = iVb ( ^ 
\7b 



iV(7) = iVb ( -L ) e- 



Cm 7min < 7 < 7b 

7b 5? 7 



The spectral indices are pi = 1.5 and p2 = 2.5. The break is at 
7b = 2 X 10", the high-energy cut-off at 7niax = 2 x 10^. The 
number density of this 'background' population is ric — 4cm~'^. 
Their total energy content is 2.2 x 10^® ergs. 

By the time when the new flare begins, i.e. the shock begins 
to cross the region and inject electrons, this pre-existing population 
has cooled to a 7max of a few times 10*, yielding a synchrotron 
peak at around 50 eV. In the observer's frame, the cooling time- 
scale for the peak of the 'background' component is of the order of 
1 day and we could think of it as due to the aging of the electron 
spectrum from a previous active phase occurred a few days earlier. 
In most recent long observin g campaigns Mrk 421 exhibited flares 



in most recent long observin g campaigns Mrk 421 e 
on about this time-scale (e.g iTakahashi et al J 20001) . 



The injection of electrons begins at t'g 



start, inj 



5 X 10^ s, with 



a power law distribution (i ]2.6.1| >. The parameters of the injected 
spectrum are: p — 1.5, 7miii = 50, 7max ~ 1.9 x 10^, I/(nj = 
5.5 X 10"" erg/s. 

The emitted -beamed- photons are integrated over the angle 
of 0.99944 < cos(6') < 0.99964, which corresponds to a Doppler 
factor of 27 < 5 < 42. 



4.4.1 Results 

In Fig. |4^-d we show a summary of the main comparisons with 
observations, as SEDs, light curves and flux-flux correlation. The 
broadband SEDs at 3 different times are shown in Fig.|4^, with X- 
ray and TeV 7-ray spectra for 2001 March 19 and historical multi- 
wavelength (from radio to TeV) data points. Corresponding SEDs 
zoomed around the X-ray band are shown in Fig. |4j). Light curves 
for 5 relevant energy bands are plotted in Fig.]?]:, while the fluxes 
in the X-ray and TeV 7-ray bands are plotted against each other 
in Fig. |4ji. About the light curves, it is important to note that the 
evolution during the first 15 kilo-seconds (— istart.inj/^) of s™" 
ulation (highlighted with grey shading) simply reflects the initial 
setup of the pre-existing background electron population, reaching 
its (approximately) steady radiative state as the blob fills with the 
radiation from all the zones, and radiation begins to escape. 

I n Figures [5^,b we show the Discrete Correlation Function 
(DCp lEdelson & KroIik|[l988 b computed between the light curves 
in two X-ray bands (2—4 and 9— 15keV) and X-ray and TeV 
(2— 4keV and > 1 TeV). They are shown for illustrative purposes, 
and no extensive statistical analysis has been performed to assess 
the uncertainty on the lag value. 

In the framework of the observational issues outlined in Sec- 
tion|4l we note that: 
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Figure 4. Summary of results of the tirst case, with pre-existing background electron population. All quoted times are in the observer's frame, (a) Broad 
band SEDs for three representative times (see labels) during the simulation of the flare. Observed spectra for high and low states during the 2001 March 19 
flare in X-ray and TeV 7-rays are plotted in cyan and blue. Other historical data points in grey or black are the same as in jpossati et al.l 12008). The black 
empty squares in the GeV 7-ray band are the FenjiiVLAT 1-year averages in five bands, (b) Zoom centered on the X-ray band, (c) Light curves at 5 different 
energy bands, on 700 seconds bins. Light curves are normalized to their peak values. In the bottom panel we show two X-ray energy bands, and in blue the 
RossiXTE/PCA 2—4 keV data. The top panel comprises the simulated TeV band light curve, to be compared with the blue data points, as well as light curves 
in optical and a band representative of the Fenni/LAI bandpass. The grey shaded area is intended to highlight the initial section of the light curve which is not 
meaningful because it corresponds to the interval during which the pre-existing electron population is being 'prepared' . The long dashed vertical grey lines 
mark the injection period. The dotted red line marks the time coiresponding to the largest cross-section of the active region along planes of equal observed 
times (see text, ^ 14.7) . The colored short thick segments mark the times corresponding to the SEDs plotted in (a) and (b). (d) The flux vs. flux plot for X-ray 
and 7-rays. Simulations data have been smoothed with a one-hour width boxcar filter Colors highlight different 10 ks time intervals. As labeled, red starts at 
t=30ks. For the red and blue lines we used a thicker trait to highlight the central time interval of the outburst. The magenta points connected by a dotted line 
show the observed coirelation of TeV flux in Crab units vs. X-ray in count rate for the 2001 March 19 flare. Because of the different units they are plotted at 
an arbitrary position. 



(i) The flare light curves are approximately symmetric for both 
the X-ray bands as well as for the TeV 7-rays. The GeV 7-ray light 
curve asymmetry reflects the relative duration of the light crossing 
times and of the cooling time-scales of the electrons emitting the 
seed photons and doing the IC scattering. For photons emitted by 
electrons (and seed photons) varying on a time-scale shorter than 
the geometric one this latter dominates the flare profile, hence mak- 



ing it symmetric. For bands whose emission processes are charac- 
terized by physical time-scales longer than the geometric ones, the 
-slower- cooling decay profile emerges. 

(ii) The amplitude correlation between X-ray and 7-ray fluxes 
is a quadratic, i.e. F-^/F^fi ~ {Fx/ Fx,o)^ with rj ~ 2, during the 
rise of the flare. Shortly after the peak the trend flattens, becoming 
linear. At this point the March 19 light curves were still showing a 
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Figure 5. Discrete cross coiTelations analysis. Our simulated cases run left to right, (a, b) for the 'background' case, (c, d) for the 'foreground' case, (e, f) for 
the 'better TeV spectrum' case. Top panels show the intra-band X-ray correlation between X-ray (2 — 4keV) and X-ray (9 — 15keV), and bottom panels the 
inter-band correlation between soft(er) X-ray (2 — 4keV) and TeV 7-rays. 



quadratic correlation, which lasted until the end of the TeV (Whip- 
ple) observational coverage (see magenta points in Fig.|4l(. 

(iii) A soft lag is clearly discernible between the different X-ray 
bands, while a similarly short hard lag is present between the 7-ray 
and the softer X-ray band (see also Figures[5^,b). In th e 2001 March 
19 fla re a short hard lag was observed in both cases jFossati et alj 
l2008h . We will discuss a possible important factor responsible for 
the soft intra-band X-ray lags and the role played by geometry and 
LCTE later, in SectionKTl 

(iv) Since the active region was previously filled by a population 
of electrons emitting a lower luminosity slowly varying SED this 
scenario easily accounts for the modest variability in the optical 
band. 

(v) While matching the observed X-ray spectra is relatively 
easy, for the 7-ray spectrum we encountered the usual challenge: 
the obse rved 7-ray spectrum is harder than what pred icted by sim- 
ulations jFossati et alj2000bl : lBlazeiowski et al.ll2005l) . 

In order to investigate these points in more details, we ex- 
plored two alternative scenarios, which we discuss in turn below. 



4.5 Case 2: injection in empty blob, with emission diluted by 
a separate steady component (foreground) 

With a similar setup we tested a scenario in which there is no back- 
ground electron population pre-existing in the blob. The steady 
broader band emission observed in the optical band is attributed to 
a component from a different region in the jet, which we will call 
'foreground' component. We assume that there is no interaction be- 
tween the two regions. The 'foreground' component is combined a 
posteriori with the radiation from the flaring blob, simply by adding 
it as a steady SED to the emission from the time dependent simula- 
tion. A very important difference with respect to the previous case 



is that photons from this component do not contribute to the IC 
emission by the freshly injected electron population. 

The volume size, geometry and Doppler factor of the active 
blob are the same as for the previous case. Because of the lack 
of extra local seed photons for the IC emission, in order to match 
the SED, in particular to boost the IC component with respect to 
the synchrotron one, it is necessary to decrease the magnetic field 
strength. The injection of electrons begins at fstart.inj = 5 x 10^ s. 
The injected distribution has a spectrum with p = 1.5, 7min ~ 50, 
7max = 1.9 X 10^ L;„j = 6 X 10^" erg/s. 

The foreground component is simulated with the same code, 
run separately. For convenience, its electrons are assumed to be in 
similar geometric and magnetic environment to the active region. 
They have a broken power law distribution, with spectral indices 
Pi = 1.5 and p2 — 2.5. The break is at 7b = 10*, the high-energy 
cut-off at 7max ~ 10^. The electron density is ric = 6cm^'^ (total 
energy content is 2.3 x 10*" ergs). These parameters for the putative 
'foreground' emission are such that its time evolution is modest on 
the time-scales in which we are interested here. 



4.5.1 Results 

The resulting SEDs, light curves, and the X-ray vs. TeV flux-flux 
correlation are shown in Fig.|6^-d, the DCFs in Fig.|5j:,d. This sce- 
nario does not reproduce the main features of the reference obser- 
vations better than the first one. 

(i) The flare is asymmetric in TeV 7-rays and in the softer X-ray 
band. It remains symmetric for harder X-rays. 

The relative length of cooling and geometric time-scales is again 
an important factor, cleanly shown by the X-ray bands. In TeV 7- 
rays however this is compounded by the effect of LCTEs. The first 
steeper rise (up to f = 30 ks) of TeV flux is driven by the increase of 
electrons as they are injected in the blob by the moving shock com- 
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Figure 6. Summary of results for the second case (empty blob, with foreground emission). All panels, colors, symbols are the same used in Fig.|4] with the 
addition in a) and b) of a cyan SED representing the 'foreground' component. In this case there is no preliminary phase to prepare the active region. 



bined with the fact that we see a larger and larger fraction of the 
blob volume, modulated by external LCTE (see ij |4.7t . This is also 
signaled by the fact that the knee occurs at around the time when 
the observer would see the largest section of the blob (red dotted 
line in Fig.|6j;). The slow rising, flat-top, phase (30 — 45ks) of the 
TeV light curve is due to the increase of seed photons available at 
each location within the blob due to diffusion from the rest of the 
blob, delayed by internal LCTE. It's a slow rise also because the 
high energy electrons responsible for most of the IC scattering to 
the TeV band are already cooling rapidly. At some point the radia- 
tion energy density in each location in the blob will stop increasing 
because enough time has passed for photons to diffuse throughout 
the blob. After that time the evolution is simply determined by par- 
ticle cooling and external LCTE. Because the electrons emitting the 
bulk of the observed TeV flux have a cooling time larger than R/c, 
in this case the TeV flare decay shape is determined by cooling 
rather than LCTE. 



(ii) As in case 1, the flux-flux amplitude correlation is repro- 
duced only partially. The trend is almost quadratic during the rising 
phase of the flare, and it turns to sub-linear on the decaying phase, 
after a short horizontal shift corresponding to the flat top of the TeV 
light curve. 

(iii) The path of the flux-flux diagram signals the presence of a 
time lag between the soft X-ray and the TeV 7-ray emission, which 
is shown in the DCF (Fig.[5jl). The TeV 7-ray lags the 2-4 keV soft 
X-ray by about 2 kilo-seconds, comparable to the observation of the 
March 19 flare as in the first scenario. Also similar to case 1 is the 
soft lag between the two X-ray bands, opposite to what observed 
on 2001 March 19. 

(iv) For what concerns the optical band, since we designed also 
this second scenario to address directly its minimal variability, it is 
not surprising that the light curve exhibits only a modest variation. 

(v) Finally, also in this scenario we have not been able to pro- 
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Figure 7. Summary of results for the case 3 (blob with pre-existing background electron population, with parameters adjusted to better match the TeV 
spectrum.) All panels, colors, symbols are the same used in Fig.|4] 



duce TeV 7-ray spectra as hard as the observations and in the end 
we limited ourselves to matching the flux level at around 1 TeV. 



4.6 Case 3: with pre-existing electron population, adjusted to 
better match the TeV spectrum 



Some of the differences with respect to the first case are ulti- 
mately due to the weaker magnetic field making synchrotron cool- 
ing time-scales ~ 60% longer: for the highest energy electrons 
emitting in X-ray and TeV r^ooi becomes longer than Rjc. It is 
worth emphasizing that the decrease of B is dictated by observa- 
tional constraints, namely the relative luminosity of the synchrotron 
and IC components and the need to compensate for the absence of 
the additional source of seed photons for IC scattering provided in 
case 1 by the co-located 'background' component. This is in fact a 
good example of how the model is globally constrained. 



As we pointed out, in the previous two cases, the simulated SED 
in the TeV 7-ray range is softer than the observed spectra. To 
try to improve the match of the TeV part of the SED, we con- 
sidered a modified version of the first scenario. We increased the 
Lorentz factor (F ~ 46) and decreased the magnetic field strength 
{B — 0.035 G), the goal being to move the inverse Compton 
peak to higher energy while leaving the synchrotron peak approx- 
imately unchanged. The parameters are: R = 1.5 x 10^^ cm, 
Z ^ 2x 10^^ cm, B = 0.035 G, T = 46. At f = the 'back- 
ground' electrons have the same broken power law distribution as 
in the first case, but with lower number density, ric = 1.56 cm"'^. 
The volume is slightly larger yielding a total energy content of 
' ergs. 

The injected electrons have a power law distribution with p = 
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1.9 X 10^ L;„j = 3.2 X W^erg/s. 
8 X 10^ s. 



1.5, 7min = 50, 7max 

Injections starts in this case at fgtart.inj 

Results are integrated over 0.99971 < cos(6') < 0.99981 
which corresponds to a Doppler factor range 37 < 5 < 57. 



4.6.1 Results 

SEDs, light curves, and X-ray vs. TeV flux-flux correlation are 
shown in Fig. [T] and DCFs in Fig. |5^,f. The itemized summary 
of the main reference observations does not show improvements 
beyond the slightly higher VHE SED peak. 

(i) Because of the larger Doppler factor the electron emitting at 
the SED peaks have lower energy, which combined with a weaker 
magnetic field results in longer cooling time-scales (see eq.ll2t. in 
turn exceeding the source crossing time. This has the effect of in- 
creasing the asymmetry of the light curves in bands whose emission 
involves lower energy electrons and/or photons. The soft X-ray and 
TeV 7-ray light curves indeed have a slowly decaying tail. 

(ii) Once again during the rising phase of the flare the 7-ray- 
X-ray correlation is approximately quadratic, until the peak of the 
X-ray light curve. After the TeV flare peak the correlation is ap- 
proximately linear, as expected when the variation in both bands 
is driven only by the cooling of the (same) electrons, because the 
IC seed photons are emitted by particles with longer cooling time- 
scale. 

(iii) The results concerning time-lags are equivalent to those of 
the other scenarios, perhaps with a hint of a smaller lag between 
TeV and softer X-ray. More extensive analysis would be necessary 
to quantify this possibility. 

(iv) The slight shift of the IC peak to higher energy enables a 
better match with the observed spectra, although the actual spec- 
tral indices of the simulated SEDs remain softer than the observed 
values. 

Further increases of the Doppler factor can still produce good 
SEDs, as long as we concurrently reduce the size of the volume. 
However, the light crossing time would rapidly become smaller 
than the observed flare duration, and it would have minimal im- 
pact on the observed phenomenology. Therefore the observed flare 
shapes must represent the true acceleration and cooling of the elec- 
trons, and the symmetry of the light curves must be caused by sim- 
ilar heating (or injection) and cooling time-scales. 



4.7 Geometric Effects on Light Curves 

There are complex geometry-related effects that have an impact 
not only on the shape of the observed light curve (e.g. its sym- 
metry), but can also leave an imprint on other observables such 
as time lags and energy-dependent flare shape. Depending on how 
the particle injection and acceleration processes are distributed spa- 
tially, differences in physical time-scales for particles of differ- 
ent energy effectively may add a further geometric effect by in- 
ducing inhomogeneities (e.g. stratification) in the source (see also 
IChiaberge & Ghisellinilll999l;ISokolov et alj2004l) . 

We would like to illustrate with an extremely simple toy- 
model some aspects of the role of the geometry of the emitting 
region, and its interplay with some of the intrinsic physical time- 
scales, responsible for the fact that the peaks of the simulated X-ray 
light curves did not correspond to either the time when the shock 
exits the active region and injection is not present anywhere any- 
more, or to the time corresponding to the largest cross-section of 
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Figure 8. Diagram illustrating the geometry of the toy-model (see text, 
ii l4.7t . The black rectangle represents a side 2D view of the cylinder, with 
R = 1 and Z = 4/3 (same aspect ratio of our simulations). The grey 
dashed line represents the 'shock front' traveUng at the speed of light 
through the cylinder, along its axis of symmetry, in the direction marked 
by the arrows. The green arrow represents the direction to the observer, 
while the green lines (solid and dotted) are planes peipendicular to the line 
of sight. In this example the viewing angle is = 70° . The red lines rep- 
resent the loci of points whose photons reach the observer simultaneously, 
taking into account the 'shock' travel distance/time until of their activation 
and the light travel time to the observer since that moment. The blue dashed 
lines illustrate this by showing three paths of equal length from the begin- 
ning of the 'flare'. The red loci form an angle 90° ~ ip/2 with the line of 
sight or, equivalently, an angle with the front face of the volume. 



the cylindrical volume along planes of equal observed times. In 
Figures|4j;,|6j;,|7};, these two times are marked as the second dashed 
grey line and the red-dotted line, respectively. Moreover, the shift 
changes with the light curve energy band as noticeable in the case 
of X-ray light curves. 

To illustrate how time shifts are caused by the different size 
of the observable regions filled with electrons contributing most of 
the emission at those frequencies, we consider a purely geometri- 
cal model solely based on the 'appearance' of slices of different 
thickness through a cylinder. 

Like in our real blob model, a shock is traveling along the axis 
of symmetry of the cylinder turning 'on' a thin local slice. Each 
point of this slice stays 'on' for a limited time. Ton. We do not con- 
sider a variation of brightness with time, just an on/off state. We 
build light curves where 'flux' is simply the size of the volume that 
is seen 'on' by the observer at any given time. The 'on' volume 
visible at each time from the observer point of view is computed 
taking into account light travel times and it cuts through the cylin- 
der along planes yielding constant arrival time to the observer. For 
an observer viewing the cylinder at an angle ip with respect to its 
axis, the loci of points whose photons he sees simultaneously are 
planes with an inclination ip/2 with respect to the face of the cylin- 
der (90° — v'/S with respect to the cylinder axis). Figure[8]shows a 
2D schematic of the geometry of the problem. 

If the cylinder was moving with Lorentz factor F, because of 
relativistic aberration at a viewing angle 6* ~ l/F we would be ob- 
serving the radiation that in the comoving frame leaves the cylinder 
'sideways', at 90° from its axis. The observed frequencies would be 
blue-shifted and times compressed, but that would be simply a scal- 
ing factor applied uniformly to them and for convenience we can 
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Figure 9. Light curves for a purely geometrical toy-model of the observations of a cyhnder seen at an angle i/p = 90° with respect to its axis. The meaning 
of the vertical lines is the same as in Fig. |4j;,|5];, [T];, except that since the 'injection' begins at t = 0, there is no first grey dashed line, (a) Three cases 
illustrating the time shift caused by the simple variation of Ton. The (light) curves are computed for Ton = 10~^ (orange), 0.5 (red), 1.0 (blue), in units of 
R/c. The orange case, with a Ton value yielding a negligible slice thickness, peaks at the expected time, (b) Attempt at matching light curves and simulation 
data without any baseline component contribution. The red curve is computed for Ton = 0.35, the blue for Ton = 0.65. The symbols are the simulated 
light curves of case 1, shown in Fig.|4j;: red empty squares for 9— 15 keV and blue triangles for 2— 4keV. (c) Same flare light curves combined with a slowly 
decaying basehne, with time-scale proportional to Ton- 



chose to use observed frequencies and times. Therefore observing 
the toy-model at^p — 90° is equivalent to observing the relativisti- 
cally moving blob at S = 1 /F, and in turn this purely geometrical 
analysis captures some of the features of the realistic model studied 
in this paper. 

The results are illustrated in Figure [5] In the first panel, we 
show the general features of the light curves obtained with this 
model, most importantly the effect of the change in the duration 
Ton- In Figure [9^ we plot the curves for three cases, showing the 
shift of the flare peak to later time as Ton increases. For very short 
Ton the maximum is reached at the expected time, that is when size 
of the plane that is the locus of points simultaneously seen by the 
observer (red lines in Fig. [Sj is the largest possible for the given 
viewing angle. In general, however, the light curve peak will be 
shifted by At = Ton/2. There is also a widening of the light curve, 
though much less noticeable than the peak shift. It is worth noting 
that at this extreme level of simplification geometric effects can not 
produce any asymmetry in the light curves. 

We tried to reproduce with this toy-model the X-ray light 
curves from the simulations of the first scenario presented here. 
The results are shown in Figure |9j),c. The central panel shows the 
curves obtained by adjusting Ton to match the peak time of the two 
X-ray light curves (plotted with symbols), yielding a very satis- 
factory result for values around 0.35 and 0.65 R/c. However the 
overall shape of the flares can not be matched well without adding 
a baseline component, to mimic in some way the fact that in the 
simulations the blob was already 'on' at a low level prior to the 
beginning of the injection. We therefore added to the toy-model 
light curves a slowly decaying component, with an initial level such 
that the combination of the two components would match the data. 
Since in the simulations the pre-existing component is left to evolve 
(cool) starting at the beginning of the active phase, here we let the 
baseline component decay too. The visual matching is not critically 
sensitive to the exact values of the decay slopes, and to make the 
model more constrained we forced them to be in a fixed ratio with 
respect to the chosen Ton, by considering that the cooling times of 
the electrons emitting the baseline photons are also related to their 



energies. For the results shown in Figure [9}; the slope is equal to 
five times Ton for both cases. 

The synchrotron cooling time-scale for electrons emitting in 
the 2 — 4keV and 9 — 15keV bands, following the approximate 
expression of eq. [12] are Tsoft — 0.73 — 1.04 R/c, Thard — 
0.38 — 0.5 R/c. Given the steepness of the X-ray spectrum the 
emission in each band is dominated by the lower energy electrons, 
hence the longer t is probably a more appropriate estimate. On 
the other hand, the above cooling time-scales only consider syn- 
chrotron cooling. Including some additional loss due to IC would 
decrease the value of t. In any case, the similarity between these 
crude estimates of cooling time-scales and the values of Ton cor- 
roborates the success of the geometrical toy-model at fitting the 
simulation light curves. 

The ability of the purely geometrical toy-model to reproduce 
the two X-ray light curves is indeed remarkable. For X-rays this is 
facilitated by the fact that the synchrotron emission is independent 
on the internal delays due to photon diffusion that affect the evolu- 
tion of the IC emission from the blob. It is not possible to apply a 
similar toy-model to the 7-ray light curves. 

It is worth noting that although this test shows how dominant 
the effect of the geometry can be in shaping the light curve, at the 
same time we need to highlight that some geometry parameters, 
such as the 'thickness' of the visible slices, are in effect determined 
by the physical conditions of the emission region. 

In this respect it is interesting to note that, at least in the 
setup of the scenarios presented in this work, despite the appar- 
ent dominance of the source geometry the effect of the energy de- 
pendent physics-induced geometrical factors is detectable. Hence 
multiwavelength datasets and time-resolved spectroscopy have the 
potential to disentangle them from the source geometry. 



5 DISCUSSION 

We introduced a coupled Fokker-Planck and Monte Carlo code al- 
lowing us to study blazar phenomenology in unprecedented detail 
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with a time dependent and multi-zone model properly taking into 
account all light travel time effects. 

We presented three test scenarios, aimed at modeling the vari- 
ability exhibited by Mrk 421 during the 2001 March 19 flare, and 
based on a relatively standard choice of parameters. The results of 
these tests are summarized in Table[2l side by side wit h the features 
observed in the actual multiwavelength observations (Fossati et alj 
Hop 8). There are a few fundamental issues that we wanted to ad- 
dress, which are common throughout the phenomenology of all 
well studied blue blazars. 

(i) The shape of the flares, often quasi-symmetric for a wide 
range of observational bands where the intensity variations are 
large (the main ones being X-ray and 7-ray). 

(ii) The characteristics of the correlation between X-ray and 7- 
ray fluxes. There has been great interest in the slope of their rela- 
tionship, in particular because of the observation of a quadratic, or 
higher order, relationship holding throughout some well sampled 
outbursts, challenging our understanding of the physical conditions 
and causes of the variability. 

(iii) The phase of the correlation between variations in different 
bands, namely the existence of time lags and their duration. 

None of the three test scenarios was able to reproduce all the char- 
acteristics of the 2001 March 19 flare. Two features have been par- 
ticularly challenging to match: the relationship between X-ray and 
TeV 7-ray fluxes on the decay phase of the flare, and the intra- 
band X-ray time lag. Moreover, the shape (symmetry) of the flare 
light curves could be reproduced only by one of the three scenarios, 
case 1. 

These aspects of phenomenology are among those more af- 
fected by the spatial extent and geometry of the source, whose in- 
fluence varies with observed energy band because of the relative 
importance of geometrical and physical time-scales. The impact of 
the geometrical factor, both due to the source intrinsic structure and 
to the stratification of properties due to the physical processes, em- 
phasizes the necessity of a code like the one we introduce here for 
modeling the variable high energy emission from blazar jets. 

The difficulty of producing a quadratic relationship between 
the fluxes in the X-ray and 7-ray during the declining phase of 
the flare may indicate that radiative cooling cannot fully explain 
the electron cooling mechanism. The delayed evolution of the seed 
photon field due to internal LCTE compounds the problem. One 
alternative possibility could be a process causing energy loss over 
a wide range of electrons energies (such that the IC seed photons 
are also affected) on very similar/same time-scale, such as adiabatic 
cooling, which could be associated with expansion of the blob, or 
particle escape. They are often invoked in qualitative discussions 
and in the context of simpler models, treated by means of some 
phenomenological prescription. The addition of such mechanisms 
to the code in a proper astrophysical way is not immediate, but we 
are working on its implementation. The escape term present in the 
Fokker-Planck equation is actually neglected for these set of simu- 
lations. In a follow up work including particle acceleration and es- 
cape, this latter seems to be eff ective and we obta in a quadratic flux 
relationship and hard lags (e.g. lChen et alj|201 ih . About the effect 
of adiabatic expansion of the emitting blob, based on the simplified 
analysis of Katarzvnski et al. (2005), Aharonian et al. (2009), ar- 
gue against its viability once the implications of this expansion on 
the magnetic field and particle cooling are taken into account. Nev- 
ertheless, its effect should be assessed with actual time-dependent 
simulations of a source of finite size. 

As to the hard intra-band X-ray lags, in Section l477l bv means 
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of a toy-model we illustrated an important factor affecting observa- 
tions of time-lags: source 'stratification' combined with LCTE in- 
duces a systematic soft lag. The simulations presented in this work 
do not include an acceleration term, hence the lack of hard lags is 
not a complete surprise. Nevertheless, the induced systematic ge- 
ometric soft lag introduces an additional constraint on viable elec- 
tron acceleration and injection scenarios. In this respect , the results 
of the above mentioned study focused on X-ray lags dChen et 
I20T H) suggest that continuous acceleration, spati ally extended (e.g. 
diffus e diffusive acceleration due to turbulence, iKatarzviiski et alj 
l2006l) . may be necessary, possibly accompanied by an achromatic 
energy loss mechanism. It is interesting that this type of scenario 
seems to be able to produce also a quadratic X-ray/TeV relation- 
ship throughout a flare. 

One of the most interesting aspects of this analysis was the 
comparison between two possible hypotheses for the presence of 
an additional component contributing to the observed SED. The 
need for two components, a flaring and a quasi-steady one, to in- 
terpret some of the observations has become more evident with the 
improvement of multiwavelength observations (for a n interesting 
decomposition of a Suzaku spectrum of Mrk 421 see lUshio et alj 
2009). Disentangling these two components is necessary in order 
to understand the nature of the transient activity whose proper- 
ties need to be seen more clearly. This decomposition might also 
yield information on the average properties of the relativistic jet, 
for which the less variable component might be more representa- 
tive. 

We considered the two simplest possibilities: i) that the sec- 
ondary component is due to a population of electrons that exists in 
the same region that will become active, and that will be affected 
by the flare and evolve with it. This pre-existing component can be 
interpreted as due to the remnants of a previous outburst, ii) That 
the secondary component is completely independent from the flare, 
and it contributes just a steady SED diluting the transient compo- 
nent from the observer point of view. Our simulations offer some 
hints as to their viability, and they favor the first type of scenario. 
One fundamental difference between the two alternatives concerns 
the production of IC emission, i.e. the TeV band. If the observed 
SED consists of the sum of two independent contributions, emitted 
by electrons at two different locations then the only seed photons 
for IC scattering will be those produced by the injected electrons 
themselves. Starting from an empty blob, the energy density of syn- 
chrotron seed photons needs some time to build up, which naturally 
results in a delay in the variation of the IC scattered 7-rays. This 
delay is caused by internal LCTE and it turned out to be quite sig- 
nificant as illustrated by the TeV light curve of case 2 (Fig. |6j;). 
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yielding a flat-top flare not seen in observations. That TeV light 
curve is in fact an excellent example of LCTEs at work and of the 
importance of a more advanced modeling code. 

Naturally, the cases presented in this paper only represent an 
initial study aimed at investigating the importance of LCTE which 
for the first time could be fully accounted for. These results can not 
be considered conclusive. Nonetheless, despite their limited scope 
they make a strong case for a true time-dependent and multi-zone 
modeling. 

The three scenarios discussed can be generally regarded as ho- 
mogeneous blob scenarios. The magnetic field is the same through- 
out the simulation volume, and isotropic. For what concerns elec- 
trons, the initial setup is homogeneous and the injection is identi- 
cal in all zones. It is, however, worth emphasizing that during the 
evolution of the simulation electrons properties become inhomoge- 
neous because of the different radiative cooling they experience in 
different zones. 

These simulations represent a first order implementation of 
a class of scenarios for blazar flares often discussed in the litera- 
ture, envisaging a shock acting on a discrete blob/shell within the 
jet. We adopted a volume with relatively symmetric aspect ratio, to 
not depart too much from the sphere 'implied' by one-zone models 
while making it possible to appreciate the effects of geometry, and 
attributed the flare to the injection of a fresh electron population. 

As discussed in Sections 14.11 14.21 and 14.31 the main physical 
parameters are fairly well constrained and the results can be re- 
garded as meaningful for what concerns the time- varying compo- 
nents of the model, as well as the nature of the secondary emis- 
sion component. The simulations presented in this paper suggest 
that a simple injection in the radiating region of particles with a 
formed spectrum produced in a separate acceleration region whose 
emission is not significant, does not provide a satisfactory match 
with some basic observational facts. Moreover, the comparison of 
the background and foreground component scenarios, in particular 
with respect to the TeV 7-ray evolution, clearly favoring the first 
one (case 1), suggests that if a flare is caused by a change affect- 
ing the electron population it may be necessary for it to happen on 
a relatively hot blob, acting on the same particles, re-accelerating 
them. This in turn would support a scenario where flares are not 
fully independent of each other but rather occur in the same region. 



5.1 Outlook 

We are working on a few directions of expansion of this investiga- 
tion. The main immediate focus is on the two unsolved issues of 
hard intra-band X-ray lags and X-ray-TeV flux correla tion and we 
have introduced particle acceleration and escape (e.g. IChen et alj 
l20TTh . 

Further areas of study concern inhomogeneities and different 
geometries, which can be easily studied with this code. We can 
introduce a spatial structure to the magnetic field, either static or 
changing according to some prescription, which could be motivated 
as caused by compression and amplification of the tangled field by 
a shock, and we plan to expand the code to deal with anisotropic 
magnetic field. We have started to simulate more 'exotic' geome- 
tries, e.g. elongated or flattened blobs, and the effect of an energy 
release in a small sub-region (bubble) embedded in the simulated 
volume. 

A particularly timely line of investigation concerns red blazars 
like FSRQs. Because the peak of their high energy component 
occurs in the less accessible MeV-GeV band, they have received 



much more limited attention, except at the time of EGRET, 
that however lacked the sensitivity to follow even bright sources 
throughout their full variability cycles. FermifLAT has changed 
the status quo by providing continuous coverage of several bright 
blazars, and even more importantly by being able to detect them 
during their more quiescent phases. Among the most remarkable 
examples there are 3C 454.3 and PKS 1510-089 that already pro- 
vided new clues about the nature of their 7-ray emission and the 
structure of the jet by the study of the correlated variations between 
7-rays and the ir synchro tron or thermal em ission jAbdo et al .120091 : 
iBonning et al . 2009; M arscher et al.ll201ol) . Our code allows us to 
simulate EC scenarios, by illuminating the active region with an 
external radiation field with full treatment of the relativistic aberra- 
tions, thus allowing us to model external components with different 
spectral, spatial and temporal properties. 
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